# This code reproduces Fig. 3, Panel B and associated ANAVAs #

#library(foreign)
library(optiscale)
library(lattice)

Vals <- read.table(file.choose(), header = TRUE)

# Q73_1 = Freedom, Q73_2 = Equality, Q73_3 = Econ. Sec.
# Q73_4 = Moral Trad., Q73_5 = Law & Ord

# Coded so #1 is most important, #5 is least important

#Vals$RankFree <- 5 - Vals$Q73_1
#Vals$RankEq <- 5 - Vals$Q73_2
#Vals$RankES <- 5 - Vals$Q73_3
#Vals$RankMT <- 5 - Vals$Q73_4
#Vals$RankLO <- 5 - Vals$Q73_5

table(Vals$RankFr)
table(Vals$RankEq)
table(Vals$RankES)
table(Vals$RankMT)
table(Vals$RankLO)

#Vals$Lib <- ifelse(Vals$Q36 < 4, 1, 0)
#table(Vals$Lib)

#Vals$Con <- ifelse(Vals$Q36 > 4, 1, 0)
#table(Vals$Con)

#Vals$Mod <- ifelse(Vals$Q36 == 4, 1, 0)
#table(Vals$Mod)

Vals$Dem <- ifelse(Vals$Party < 4, 1, 0)
Vals$Rep <- ifelse(Vals$Party > 4, 1, 0)
Vals$Ind <- ifelse(Vals$Party == 4, 1, 0)

Vals$Lib <- ifelse(Vals$Ideol < 4, 1, 0)
Vals$Con <- ifelse(Vals$Ideol > 4, 1, 0)
Vals$Mod <- ifelse(Vals$Ideol == 4, 1, 0)

Vals$StrongSortedDem <- ifelse(Vals$Party == 1 & Vals$Ideol == 1, 1, 0)
table(Vals$StrongSortedDem)
Vals$WeakSortedDem <- ifelse(Vals$Dem == 1 & Vals$Lib == 1 & Vals$StrongSortedDem == 0, 1, 0)
table(Vals$WeakSortedDem)
Vals$UnsortedDem <- ifelse(Vals$Dem == 1 & Vals$Lib == 0, 1, 0)
table(Vals$UnsortedDem)
Vals$StrongSortedRep <- ifelse(Vals$Party == 7 & Vals$Ideol == 7, 1, 0)
table(Vals$StrongSortedRep)
Vals$WeakSortedRep <- ifelse(Vals$Rep == 1 & Vals$Con == 1 & Vals$StrongSortedRep == 0, 1, 0)
table(Vals$WeakSortedRep)
Vals$UnsortedRep <- ifelse(Vals$Rep == 1 & Vals$Con == 0, 1, 0)
table(Vals$UnsortedRep)

Vals$SortedParty <- NA
Vals$SortedParty[Vals$StrongSortedDem == 1] <- "SSD"
Vals$SortedParty[Vals$WeakSortedDem == 1] <- "WSD"
Vals$SortedParty[Vals$UnsortedDem == 1] <- "UD"
Vals$SortedParty[Vals$Ind == 1] <- "I"
Vals$SortedParty[Vals$StrongSortedRep == 1] <- "SSR"
Vals$SortedParty[Vals$WeakSortedRep == 1] <- "WSR"
Vals$SortedParty[Vals$UnsortedRep == 1] <- "UR"
table(Vals$SortedParty)

###
# Missing Data Stuff
###

Vals$Missing1 <- NULL

Vals$CaseID <- seq.int(nrow(Vals))

Vars <- c("SortedParty", "RankEq", "RankLO", "RankFr", "RankES", "RankMT", "CaseID")

FullData <- Vals[Vars]
FullData <- na.omit(FullData)

FullData$AllSame  <- ifelse(FullData$RankFr == FullData$RankEq &
                              FullData$RankFr == FullData$RankES &
                              FullData$RankFr == FullData$RankLO &
                              FullData$RankFr == FullData$RankMT, 1, 0)
table(FullData$AllSame)

FullData <- FullData[which(FullData$AllSame == 0), ]

Vals$Missing <- Vals$CaseID %in% FullData$CaseID
Vals$Missing1 <- ifelse(Vals$Missing == "FALSE", 1, 0)
table(Vals$Missing1)

Vals$Female <- NA
Vals$Female[Vals$Q46 == 1] <- 0
Vals$Female[Vals$Q46 == 2] <- 1

Vals$Feminine <- Vals$Q68

Vals$Income <- Vals$Q49

Vals$Educ <- Vals$Q50

Vals$Educ2 <- NA
Vals$Educ2[Vals$Q50 == 1] <- 1 # Less than HS
Vals$Educ2[Vals$Q50 == 2] <- 2 # HS Diploma
Vals$Educ2[Vals$Q50 == 3] <- 3 # Some College
Vals$Educ2[Vals$Q50 == 4] <- 3 # 2-year Degree
Vals$Educ2[Vals$Q50 == 5] <- 4 # 4-year Degree
Vals$Educ2[Vals$Q50 == 6] <- 5 # Beyond 4-year degree
Vals$Educ2[Vals$Q50 == 7] <- 5 # Beyond 4-year degree

Vals$Relig <- 6 - Vals$Q52

Vals$Age <- Vals$Q87

Vals$Nonwhite <- ifelse(Vals$Q47 == 1, 0, 1)
table(Vals$Nonwhite)



FullData$AllSame <- NULL

ValsOnly <- c("RankEq", "RankLO", "RankFr", "RankES", "RankMT")
ValsData <- FullData[ValsOnly]


ValsData[1:10, ]



###

###
###   Standardize data within rows
###

vals.std <- t(apply(ValsData, 1, scale))

vals.std[1:10,]

###
###   Initialize the matrix of optimally-scaled
###   values, the previous fit, the iteration
###   number, and the improvement in fit over
###   the previous iteration
###

vals.os <- vals.std

prev.fit2 <- 0

niter <- 0

improve = 1

###
###   Start iterations
###

while (improve > .01 & niter <= 25) {
  niter <- niter + 1
  
  ###
  ###   Perform SVD on OS version of rank-orders
  ###
  
  decomp <- svd(vals.os)
  
  ###
  ###   Calculate 2-dimensional goodness-of-fit
  ###   and improvement in fit over previous iteration
  ###
  
  d.sqd <- decomp$d^2
  
  fit2 <- sum(d.sqd[1:2]) / sum(d.sqd)
  
  fitvector <- cumsum(d.sqd) / rep(sum(d.sqd), length(d.sqd))
  
  improve <- fit2 - prev.fit2
  
  ###
  ###   Create iteration history
  ###
  
  if (niter == 1) {
    history <- c(niter, fit2, improve)
  }
  if (niter > 1)  {
    history <- rbind(history, c(niter, fit2, improve))
  }
  
  ###
  ###   Obtain terminal points of vectors for row objects
  ###   (respondents in this case), calculate predicted
  ###   ranks, norm the row object vectors to unit length
  ###
  
  subj.c1 <- decomp$u[, 1:2] %*% diag(decomp$d[1:2])
  
  pred.vals <- subj.c1 %*% t(decomp$v[,1:2])
  
  sumsqd <- apply(subj.c1^2, 1, sum)
  
  root.sumsqd <- (matrix(sumsqd, nrow = length(sumsqd), ncol = 1)) ^ .5
  
  subj.coord <- subj.c1 / (root.sumsqd %*% matrix(1, nrow = 1, ncol = 2))
  
  ###
  ###   Obtain new optimally scaled data values,
  ###   using predicted ranks from fitted model
  ###
  
  for (i in 1:nrow(vals.os)) {
    opped <- opscale(x.qual = vals.std[i,],
                     x.quant = pred.vals[i,],
                     level = 2,
                     process = 1,
                     rescale = T)
    vals.os[i,] <- opped$os
  }
  
  ###
  ###   Set current fit to previous fit for next iteration
  ###
  
  
  prev.fit2 <- fit2
  
}

###
###   Print iteration history
###

history

###
###   The "fitvector" object shows the R-squared
###   for the solution in each dimensionality
###

fitvector

###
###   Plot value points
###

val.coords <- as.data.frame(decomp$v[, 1:2])

rownames(val.coords) <- colnames(ValsData)

val.coords

xyplot(V2 ~ V1, val.coords,
       aspect = 1,
       xlim = c(-.99, .99),
       ylim = c(-.99, .99),
       panel = function (x, y) {
         panel.xyplot(x, y, col = "black")
         panel.text(x, y, labels = rownames(val.coords),
                    pos = 1, cex = .75)
       }
)

###
###   Plot terminal points of row object vectors
###   in same space as value points
###

subj.coord <- as.data.frame(subj.coord)

set.seed(123)
xyplot(V2 ~ V1, val.coords,
       aspect = 1,
       xlim = c(-1.2, 1.2),
       ylim = c(-1.2, 1.2),
       panel = function (x, y) {
         panel.xyplot(x, y, col = "black", pch = 16)
         panel.xyplot(jitter(subj.coord$V1, amount = .05), 
                      jitter(subj.coord$V2, amount = .05), col = "black")
         panel.text(x, y, labels = rownames(val.coords),
                    pos = 1, cex = .75)
       }
)

###
###   Calculate overall mean direction and
###   mean resultant length
###

meand1 <- mean(subj.coord$V1)

meand2 <- mean(subj.coord$V2)

meand1 

meand2

mean.length <- (meand1^2 + meand2^2)^.5

mean.length

###
###   Plot mean vector along with 
###   individual vectors and value points
###

set.seed(123)
xyplot(V2 ~ V1, val.coords,
       aspect = 1,
       xlim = c(-1.2, 1.2),
       ylim = c(-1.2, 1.2),
       panel = function (x, y) {
         panel.xyplot(x, y, col = "black", pch = 16)
         panel.xyplot(jitter(subj.coord$V1, amount = .05), 
                      jitter(subj.coord$V2, amount = .05), col = "black")
         panel.text(x, y, labels = rownames(val.coords),
                    pos = c(4,1,3,3,2), cex = .75)
         panel.arrows(0, 0, meand1, meand2, angle = 15, length = .15)
       }
)

### Plotting Vectors ###

### merging some data ###
# val.coords is vector model results, plotting points for values in 2d space
# vals is TESS data with NaNs omitted (see VectorModel.r); n = 622
# subj.coord is subject coordinates in 2d space (see VectorModel.r)

demos <- c("SortedParty")

vals2 <- FullData[demos]

vals3 <- data.frame(cbind(vals2, subj.coord))

### Graph for Party ###
# Creating 3 Category Party Var #

#vals3$party3 <- "Independent"
#vals3$party3[vals3$party >= 1 & vals3$party <= 2] <- "Republican"
#vals3$party3[vals3$party >= 6 & vals3$party <= 7] <- "Democrat"

# Now calculate mean vectors for respondents in the three categories. Also
# calculate mean resultant lengths for the category vectors.

mean.pty.d1 <- tapply(vals3$V1, vals3$SortedParty, mean, na.rm = T)

mean.pty.d2 <- tapply(vals3$V2, vals3$SortedParty, mean, na.rm = T)

pty.mean.lengths <- (mean.pty.d1^2 + mean.pty.d2^2)^.5

# All Respondents
g2 <- xyplot(V1 ~ V2, data = vals3,
             aspect = 1, main = "All Respondents",
             panel = function (x, y) {
               panel.xyplot(val.coords$V1, val.coords$V2, pch = 4, 
                            cex = .5, col = "black")
               panel.text(val.coords$V1, val.coords$V2, 
                          labels = c("eq", "lo", "free",  
                                     "es", "mt"),
                          pos = c(2, 1, 3, 4, 1), cex = .99)
               panel.segments(rep(0, 7), rep(0, 7), 
                              mean.pty.d1[c(1,2,3,4,5,6,7)], mean.pty.d2[c(1,2,3,4,5,6,7)], lwd = 1.5)
               panel.text(mean.pty.d1[c(1,2,3,4,5,6,7)], mean.pty.d2[c(1,2,3,4,5,6,7)], 
                          labels = c("I", "SSD", "SSR", "UD", "UR", "WSD", "WSR"), pos = c(2,2,2,2,2), cex = 1.1)
             },
             xlab = "",
             ylab = "",
             scales = list(draw = FALSE)
)

###

g3 <- xyplot(V1 ~ V2, data = vals3,
             aspect = 1, 
             #main = "By Party and Sorting",
             panel = function (x, y) {
               panel.xyplot(val.coords$V1, val.coords$V2, pch = 16, 
                            cex = .5, col = "black")
               panel.text(val.coords$V1, val.coords$V2, 
                          labels = c("Eq", "LO", "Free", "ES", "MT"),
                          pos = c(2, 2, 4, 4, 4), cex = .75)
               panel.segments(rep(0, 7), rep(0, 7), 
                              mean.pty.d1[c(1,2,3,4,5,6,7)], mean.pty.d2[c(1,2,3,4,5,6,7)], lwd = 1.5)
               panel.text(mean.pty.d1[c(1,2,3,4,5,6,7)], mean.pty.d2[c(1,2,3,4,5,6,7)], 
                          labels = c("Ind", "SSD", "SSR", 
                                     "UD", "UR", "WSD", "WSR"), pos = c(4,1,3,4,4,4,4), cex = 1)
             },
             xlab = "",
             ylab = "",
             scales = list(draw = FALSE)
)

### Just Repubs

meand1 <- mean(vals3$V1[vals3$SortedParty == "SSR" | vals3$SortedParty == "UR" |
                          vals3$SortedParty == "WSR"])

meand2 <- mean(vals3$V2[vals3$SortedParty == "SSR" | vals3$SortedParty == "UR" |
                          vals3$SortedParty == "WSR"])

mean.res.length <- (meand1^2 + meand2^2)^.5

res.length <- length(vals3$V1[vals3$SortedParty == "SSR" | vals3$SortedParty == "UR" |
                                vals3$SortedParty == "WSR"]) * mean.res.length

ideol.mean.lengths <- c(.522, .312, .432)

incparty.between <- sum(ideol.mean.lengths * 
                          as.vector(table(vals3$SortedParty[vals3$SortedParty == "SSR" | 
                                                              vals3$SortedParty == "UR" |
                                                              vals3$SortedParty == "WSR"]))) - 
  res.length

incparty.within <- length(vals3$V1[vals3$SortedParty == "SSR" | 
                                     vals3$SortedParty == "UR" |
                                     vals3$SortedParty == "WSR"]) - 
  sum(ideol.mean.lengths * 
        as.vector(table(vals3$SortedParty[vals3$SortedParty == "SSR" | 
                                            vals3$SortedParty == "UR" | 
                                            vals3$SortedParty == "WSR"])))

dfnum <- 2

dfdenom <- length(vals3$V1[vals3$SortedParty == "SSR" | vals3$SortedParty == "UR" |
                             vals3$SortedParty == "WSR"]) - 3

incparty.ms.between <- incparty.between / dfnum

incparty.ms.within <- incparty.within / dfdenom

F.incparty <- incparty.ms.between / incparty.ms.within

obs.prob.incparty <- 1 - pf(F.incparty, dfnum, dfdenom)

# F on 2 and 434 df = 4.70 p = .009

### Just Dems

meand1 <- mean(vals3$V1[vals3$SortedParty == "SSD" | vals3$SortedParty == "UD" |
                          vals3$SortedParty == "WSD"])

meand2 <- mean(vals3$V2[vals3$SortedParty == "SSD" | vals3$SortedParty == "UD" |
                          vals3$SortedParty == "WSD"])

mean.res.length <- (meand1^2 + meand2^2)^.5

res.length <- length(vals3$V1[vals3$SortedParty == "SSD" | vals3$SortedParty == "UD" |
                                vals3$SortedParty == "WSD"]) * mean.res.length

ideol.mean.lengths <- c(.447, .343, .613)

incparty.between <- sum(ideol.mean.lengths * 
                          as.vector(table(vals3$SortedParty[vals3$SortedParty == "SSD" | 
                                                              vals3$SortedParty == "UD" |
                                                              vals3$SortedParty == "WSD"]))) - 
  res.length

incparty.within <- length(vals3$V1[vals3$SortedParty == "SSD" | 
                                     vals3$SortedParty == "UD" |
                                     vals3$SortedParty == "WSD"]) - 
  sum(ideol.mean.lengths * 
        as.vector(table(vals3$SortedParty[vals3$SortedParty == "SSD" | 
                                            vals3$SortedParty == "UD" | 
                                            vals3$SortedParty == "WSD"])))

dfnum <- 2

dfdenom <- length(vals3$V1[vals3$SortedParty == "SSD" | vals3$SortedParty == "UD" |
                             vals3$SortedParty == "WSD"]) - 3

incparty.ms.between <- incparty.between / dfnum

incparty.ms.within <- incparty.within / dfdenom

F.incparty <- incparty.ms.between / incparty.ms.within

obs.prob.incparty <- 1 - pf(F.incparty, dfnum, dfdenom)

# F on 2 and 561 df = 1.58, p = .21
